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INTRODUCTION 


The  U.  S.  Army  Atmospheric  Sciences  Laboratory  (ASL)  maintains  a  systems-level 
code,  IMTURB1,  which  is  used  to  calculate  propagation  statistics  for  optical  systems 
operating  within  the  atmospheric  boundary  layer.  The  code  contains  an  environment 
module  and  a  propagation  module.  The  environment  module  estimates  the  height- 
dependent  turbulence  level;  the  inner  scale,  and  the  outer  scale  as  a  function  of  me¬ 
teorological  variables.  The  propagation  module  uses  the  environmental  parameters  to 
predict  the  effects  of  propagation  along  extended  paths  within  the  turbulent  region. 
The  propagation  model  is  based  on  a  1976  propagation  code  developed  by  Fried  [1]. 
The  basic  parameters  computed  are  (1)  the  log-amplitude  variance,  (2)  the  receiver  co¬ 
herence  diameter  for  long-  and  short-term  averages,  (3)  the  isoplanatism  effective  path 
length,  and  (4)  the  scintillation  averaging  length  [2]. 

Fried’s  model  is  based  on  the  Rytov  approximation,  which  formally  restricts  the 
results  to  small  amplitude  fluctuations.  Nonetheless,  the  Rytov  approximation  is  at¬ 
tractive  because  the  propagation  parameters  can  be  computed  in  terms  of  comparatively 
simple  path  integrals.  Over  the  past  few  decades,  however,  considerable  progress  has 
been  made  in  the  iheory  of  wave  propagation  in  random  media,  particularly  under  the 
narrow-angle  scatter  approximation.  Thus,  it  is  no  longer  necessary  to  accept  the  limi¬ 
tations  of  the  Rytov  approximation.  Indeed,  the  linear  systems  approach  to  modeling 
imaging  systems  leads  to  image  degradation  measures  that  are  expressed  in  terms  of  the 
mutual  coherence  function,  which  is  amenable  to  exact  evaluation.  The  linear  systems 
model  is  reviewed  in  the  appendix  to  this  report. 

The  report  proper  describes  a  new  propagation  module  that  will  allow  the  IM- 
TURB1  code  to  accommodate  essentially  unrestricted  propagation  disturbance  levels, 
while  preserving  the  compact  path-integral  formulation.  In  work  performed  under  a 
previous  contract,  the  Shkarofsky  spectral  density  function  model  was  used  to  derive  an 
analytic  form  for  the  phase  structure  function  that  retains  explicit  inner  and  outer  scale 
cutoff  parameters  [3].  This  provides  an  efficient  method  for  computing  the  complete 
mutual  coherence  function  (MCF)  in  place  of  the  single-parameter  characterization  that 
is  currently  used  in  the  IMTURB1  model.  For  plane  and  spherical  waves  this  is  easily 
done  because  the  MCF  admits  a  simple  analytic  characterization  in  terms  of  a  path 
integral  over  the  phase  structure  function.  For  beam  waves,  the  corresponding  path 
integral  is  nested  in  a  two-dimensional  integral  that  must  be  evaluated  numerically 
or  approximated.  An  approximation  proposed  by  Ishimaru  is  used  for  the  proposed 
propagation  module. 

To  provide  a  consistent  framework  for  computing  the  effects  of  long  and  short  in¬ 
tegration  times,  a  scheme  proposed  by  Fante  has  been  implemented  in  which  the  low- 
frequency  spectral  contribution  is  removed  before  the  structure  function  is  computed. 
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Thus,  the  same  basic  theoretical/computational  framework  can  be  used  for  all  the  prop¬ 
agation  parameters.  The  short-term  results  are  incorporated  as  correction  factors  to  the 
long-term  structure  function  with  the  necessary  numerical  integrations  approximated  by 
a  set  of  polynomial  functions.  This  allows  efficient  computation  within  the  propagation 
module. 

For  completeness  we  have  repeated  some  earlier  material  describing  the  general 
solution  to  the  parabolic  wave  equation  for  the  mutual  coherence  function  (Section  11) 
and  the  Shkarofsky  spectral  form  that  is  used  to  model  the  phase  structure  function 
(Section  III).  We  then  review  the  specific  forms  of  the  mutual  coherence  function  that 
are  proposed  for  plane,  spherical,  and  beam  waves  (Section  IV).  Section  V  describes 
the  modifications  for  short  integration  times.  The  results  form  a  complete  and  consis¬ 
tently  formulated  propagation  module,  which  can  be  integrated  into  the  IMTURB  code 
framework. 

To  provide  some  guidelines  for  interpreting  the  theoretical  results,  a  numerical  simu¬ 
lation  has  been  implemented  as  was  recently  done  by  Martin  and  Flatte  [4].  To  translate 
the  results  to  imaging  systems,  however,  the  simulated  beam  wave  field  was  refocused 
to  a  point.  Several  measures  of  the  averaged  refocused  beam  degradation  were  then 
computed  to  evaluate  both  short-term  and  long-term  beam  degradations.  These  results 
are  described  in  Section  VI. 
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II  BACKGROUND 


Propagation  of  light  in  the  atmospheric  boundary  layer  is  governed  by  the  parabolic 
wave  equation 

+  v' t,'(r)  +  =  °.  (i) 

where 

u{  r)  =  U'(r)exp{ikz}.  (2) 

For  scintillation  studies,  one  typically  assumes  that  the  relative  permittivity,  e1(  consists 
of  a  locally  invariant  component  plus  a  purely  random  perturbation;  however,  it  is 
instructive  to  accommodate  an  explicit  secular  variation,  whereby 

^1  =  Ci  +  &,  (3) 

with  ei  representing  the  slowly  varying  background  component.  As  a  first  step  in 
isolating  the  effects  of  Sj,  let 


U'(t)  =  U(  r)exp{i'V>(r)}.  (4) 

One  approach  attempts  to  identify  tJj  as  a  ray-optics  component,  in  which  case  it  has 
both  real  and  imaginary  components  as  does  the  diffractive  component  U .  It  is  simpler, 
however,  to  define  if)  as  the  solution  to 


dij)(  r) 
~  dz 


By  substituting  (3)  and  (4)  into  (1)  and  using  (5),  it  is  readily  shown  that 


=  V\U  -  {V2±j>)U  =  k2SeU. 


(5) 


(6) 


Now  assume  that  \V2xi>\  «  \VlU/U\t  whereby  the  ^'-dependent  term  in  (6)  can  be 
neglected,  and  U  itself  satisfies  the  parabolic  wave  equation. 

It  follows  that  solutions  to  the  parabolic  wave  equation  based  only  on  the  homo¬ 
geneous  random  field  &  omit  the  effects  of  slow  phase  variation  induced  by  In  an 
imaging  system,  this  phase  variation  will  cause  a  small,  slowly  varying  displacement  of 
the  image.  Such  effects,  however,  are  often  attributed  to  the  low-frequency  component 
of  the  homogeneous  turbulence.  Indeed  Fried  has  used  the  Rylov  phase  to  compute 
the  mean-square  linear  slope  for  homogeneous  turbulence  as  a  measure  of  the  average 
image  displacement.  Based  on  the  results  presented  in  Section  VI,  however,  we  believe 
that  this  component  is  too  small  to  account  for  the  observed  beam  wander  for  weak 
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propagation  disturbances.  To  properly  interpret  the  beam  wander,  the  inhomogeneous 
contribution  must  be  modeled  separately  and  incorporated  explicitly  as  in  (4). 

Either  way,  the  most  important  quantity  for  predicting  the  systems  effects  of  atmo¬ 
spheric  turbulence  is  the  mutual  coherence  function,  which  is  defined  as 

=<U(pc  +  Ap/2,z)U-(p .  -  Ap/2,z)>,  (T) 

where  the  angle  brackets  denote  ensemble  averaging,  A p  =  //  p,  and  pc  —  ( p  h  p')/ 2. 
Under  conditions  that  are  readily  satisfied  for  most  applications,  it  can  In;  shown  that 
T  satisfies  the  differential  equation 

{2ikFz +  <v’.  -  vi.)  +  y(/1(0)  '  ■4(A,,))} r(A/,;  =  °'  W 

where  j4(A/?)  is  the  autocorrelation  function  of  the  refractive  index  perturbation.  Thus, 
insofar  as  the  scattering  medium  is  concerned,  it  is  only  necessary  to  model  the  structure 
function  De(Ap)  =  2(^4(0)  —  A(Ap)).  Closed-form  solutions  can  be  obtained  for  point 
sources  and  the  limiting  case  of  plane  waves.  Good  approximations  are  also  available 
for  beams  as  summarized  in  Section  IV. 2. 
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Ill  PHASE  STRUCTURE  FUNCTION 


III.l  General  Results 

As  discussed  in  Section  II,  the  solutions  to  the  parabolic  wave  equation  for  the  complex 
moments  cf  U  can  be  formulated  in  terms  of  the  phase  structure  function.  Thus,  the 
phase  structure  function  is  fundamental  in  any  propagation  modeling  efTorl.  Consider  a 
rectangular  coordinate  system  with  the  z  axis  along  the  direction  of  propagation.  The 
phase  path  integral  for  a  distance  lp  is  given  as 

&f>=k[  &(p,z)dz.  (9) 

Jo 

The  phase  autocorrelation  function  can  be  computed  from  (9)  as 

<W>=  k2  /*'  f’  <&{p,z)&(p\z')>  dzdz'.  (10) 

Jo  Jo 

For  homogeneous  statistics, 

<&&'>=  B<(Ap,Az).  (II) 

By  using  (11)  and  performing  a  straightforward  change  of  variables  in  (10),  one  obtains 

<&f>&f>'>=  k2lp  f'  (1  -  |Az|//p)  Bt{Ap,  Az)  dAz.  (12) 

For  almost  all  applications,  the  correlation  distance  along  z  is  small  compared  to  lp,  in 
which  case  the  integral  in  (12)  can  be  replaced  by 

A{Ap)=  [°°  Bc(Ap,Az)dAz.  (13) 

Let  us  now  define  the  spectral  density  function  (SDF)  as 

$(K,A:r)  =  J  J  J  Bt(Ap,  Az)  exp{-i(K  •  Ap  -|  kzAz)}  dApdz.  (14) 

The  overbar  is  used  to  denote  SDFs  normalized  as  in  (14)  rather  than  as  in  the  alter- 
native  definition  that  places  the  (27t)3  term  in  the  denominator  of  the  complementary 
spectral  integral.  If  the  SDF  is  isotropic,  it  is  readily  shown  that 

A(Ap)  =  (2tt)2  r  J»{KAp)$€(K)KdK.  (15) 

Jo 
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The  refractive  index  n  is  related  to  e  as 


n  =  y^  =  n\/T+li  +  •  (16) 

It  follows  that  4/2„  =  -R*.  Thus,  the  structure  function  for  the  refractive  index  is  given 
as 


fn(Arf 

=  i(^(0)  -  X(A,)) 

=  2n7  r(l-J0(KAp))$c{K)KdK. 

Jo 

(17) 

It  is  convenient  to  define  a 

spectral  shape  function  Q(q)  such  that 

(18) 

With 

(19) 

it  follows  that 

JJJ$n(k)dk=<8l2>. 

(20) 

We  note  that  the  structure  function  occurs  naturally  in  gaussian  random  field  theory. 
Let  U(p)  =  exp{t$(p)},  where  &f>  is  a  gaussian  process.  It  follows  from  the  properties 
of  gaussian  variables  that 


<{/(/'*>  =  exp{— ^  <(Sf>{p)  -  fy(p'))2>} 

=  cx?{-k3lpDn(AP)}.  (21) 

Indeed,  (21)  is  the  mutual  coherence  function  for  a  plane  wave  propagating  in  a  ho¬ 
mogeneous  random  medium;  however,  this  occurs  because  the  diffraction  effects  that 
modulate  the  complex  diffracted  wave  field  average  to  unity  in  the  average  coherence 
computation. 


III.2  Power-law  Spectral  Models 

The  modified  Kolmogorov  spectrum  is  defined  in  terms  of  the  structure  constant  C* 
and  the  inner  l0  and  outer  L0  scale  sizes  as 

i„(k)  =  0M3C2n(kl  +  k2)-"/‘  «cp{-*!/0.  (22) 
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where 


K  =  5.92/Zo  (23) 

kL  =  l/Io-  (24) 

The  full  wave  vector  is  conveniently  specified  in  terms  of  its  transverse  and  axial  com¬ 
ponents  as  k  =  (K,  kt).  Bold-face  symbols  indicate  vectors.  The  magnitudes  of  the 
vector  components  will  be  indicated  by  the  corresponding  ordinary  symbol.  A  problem 
with  (22)  is  that  it  does  not  admit  a  closed  form  for  A(Ap).  A  more  convenient  form 
was  developed  by  Shkarofsky  [5].  It  will  be  used  in  a  slightly  modified  form,  which  is 
summarized  in  Table  2  of  [6],  namely, 

I - 7-— — -(v+l/2)  Kv+i/2  (2 

Q(q)  =  (2tt  ot,otL)3/2J  1  -f  — (25) 

v  *  Kv-i 

where  Kv(x)  is  the  modified  Bessel  function  of  fractional  order,  it  will  also  be  convenient 
to  define 


"  =  [qQ{q)^r 

*-•  (2s) 


It  can  be  shown  that 


m  =  A(yym 

=  [  iMiy)QM  J«/  / 

Jo  Jo 

L  y1  *->/>  (2S \A+" 


(ss) 


=  4  1  + 


Thus,  Q,  R  or  and  A  have  self-similar  Fourier  transforms. 

For  small  values  of  x,  i(2x)  «  |r(i/  “  l)(at#/ax,)“*'+1.  For  all  applications  of 
interest,  a,/a^  <<  1.  Thus,  for  q  «  1/a,  (25)  simplifies  to 


8jr3'J(aI,/v^)!(1“')r(i,  +  1/2)  (  2 

r(»  - 1)  U 


GM 


—u—l/2 
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Equation  (28)  can  be  put  in  the  same  form  as  (22)  if  (28)  is  first  divided  by  (27r)3  to 
maintain  consistent  2n  conventions.  We  then  let  ai  =  y/2L0  and 


<Sn2> 


{L0)^-v)T(u  +  l/2) 

7r3/Jr(j/  -  1) 


0.033(7’, 


(29) 


or 


7rs/JC,r(y  - 1) 

r(i/  +  i/2)921'-2’ 


(30) 


where  q0  =  Lq1  and  C,  =  0.033C12.  By  noting  that  .4(0)  =  n  <&2>,  we  can  obtain  a 
consistent  expression  for  Dt(y)  or  Dn(y).  The  complete  model  is  summarized  in  Tabic  1 
in  terms  of  the  basic  parameters  (72,  lo,  Lo,  and  the  11/3  power,  which  corresponds  to 
v  =  4/3. 


Figure  1  shows  a  comparison  of  the  Kolmogorov  and  Shkarofsky  spectra'  forms  for 
a  fixed  inner  scale  of  1  cm  and  a  typical  outer  scale  range.  The  Shkarofsky  model  was 
evaluated  from  Table  1.  Standard  subroutines  were  used  to  evaluate  the  Bertel  and 
gamma  functions.  It  can  be  seen  that  significant  differences  occur  only  in  the  regime  of 
the  inner  scale  cutoff.  Because  the  detailed  mathematical  form  in  the  dissipation  regime 
is  not  known  precisely,  these  differences  are  unimportant.  Figure  2  shows  the  Shkarofsky 
form  of  the  structure  function  as  summarized  in  Table  1.  Although  the  Shkarofsky 
form  is  convenient  for  modeling,  it  can  be  unwieldy  for  analytic  computation.  Thus, 
the  approximations  are  often  used. 

In  the  inertial  subrange  Lo  <  q  <  lo,  $n  has  the  power-law  form  C,q~u ^3.  By 
carefully  evaluating  lim^^o,,  Dn(y)  it  can  be  shown  that  the  structure  function  admits 
the  complementary  power-law  form 


Dn(y) 


47raC',r(3/2  -  v) 

r(«/  +  l/2)(2i/-l)2J‘'-jl/ 


2t/-l 


(31) 


Note  that  (31)  does  not  depend  on  L0.  Carrying  out  a  Taylor  series  expansion  of  Dn(y ) 
for  small  y,  using  (31)  for  intermediate  y,  and  the  saturation  value  of  Dn(y)  for  large 
y ,  it  follows  that 


Dn{y)  = 


47r2£  I(a/2-r)ira  2 

47r2Q  _ E(3Z2zid_v2^-l 

™  °8r(y+l/2)(2i/-l)2J‘'-jy 

4it2C  -liazim 


y  <  lo 
lo  <  y  <  Lo 
y  >  L0 


(32) 
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Evaluating  the  coefficients  for  the  Kolmogorov  value  u  =  4/3,  we  obtain 


Mv) 


2.024CS4  y  <  (o 

]0 

2.913Cjy5/3  l„<y  <  i0 
1.563C3Lo/3  y  >  Lo 


(33) 


The  asymptotic  and  saturation  forms  of  (33)  are  identical  to  the  corresponding  forms 
from  (20-78)  in  Ishimaru  [7]  when  allowance  is  made  for  the  factor  of  2  difference  between 
-4(0)  —  j4(y)  and  Dn(y)  from  (17).  The  small  y  form  given  by  Ishimaru  is  about  62% 
larger,  which  reflects  the  differences  between  the  Shkarofsky  and  Kolmogorov  forms. 
Equation  (33)  is  plotted  in  Figure  3.  It  can  be  seen  that  the  quadratic  and  asymptotic 
regimes  are  narrowly  confined,  and  results  based  on  these  approximations  must  be 
restricted  accordingly. 
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Table  1.  Summary  of  Shkarofsky  Form  of  Modified  Kolmogorov  Spectrum. 


Shkarofsky  Model  u  =  4/3 


».(«)  =<*>2> 

/ - 7 — u -(-+ 1/2) If. II/! (W \A 

Q(g)  =  (2na9aL)^2sjl  +  ^-  - ^ 

-  r(x/+l/2) 


1- 

llT 


Cs  =  0.033C“ 


if  wo  fir  i  fi 

/1+4) 

yTV~l  I'  A*'-1/2  1 

V1  +  2aJ  tf„ 

<fa2>=  qo=r-‘ 


f(H-l/2K 

OLi  —  \[2Ij§  Ois  —  Iq/^/Ti 

C(y)  =<6n2>  R(y) 


ja 
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Phi (k)/Cn**2 


Kolmogorov  and  Shkarofsky  Spectra 


-3  -2  -1  •  1  2  3 

10  10  10  10  10  10  10 

k  (1/niiUr) 


Figure  1.  Comparison  of  Kolmogorov  and  Shkarofsky  spectral  models. 
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D(y)/Cn**2 


ShK&rofsky  Structure  Func t 


lo(m)  =  0.005 
2  <  Lo ( m  J  <  9 


y  (meter! 


Figure  2.  Shkarofsky  form  of  structure  function 


D(y)/Cn**2 


Kolmogorov  Structure  Function 


y  (motor) 

Figure  3.  Approximate  forms  for  structure  function  regimes. 
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IV  MUTUAL  COHERENCE  FUNCTION 


The  results  summarized  in  this  section  are  taken  from  Chapter  20-6  of  Ishimaru  [7], 
where  detailed  derivations  and  an  extensive  bibliography  can  be  found.  As  noted  in 
Section  II,  all  solutions  to  the  differential  equation  (8)  derive  their  dependence  on  the 
propagation  medium  from  the  structure  function  for  the  refractive  index  Dn(y),  where 
y  represents  the  magnitude  of  A p.  In  a  plane  normal  to  the  propagation  direction,  tue 
structure  function  is  isotropic. 

The  receiver  coherence  diameter,  r0,  which  is  defined  as  the  path  separation  in  the 
receiver  plane  over  which  the  spatial  coherence  of  the  complex  field  decreases  by  the 
factor  e-1,  is  a  useful  parameter  for  evaluating  the  effects  of  atmospheric  turbulence. 
In  the  current  version  of  IMTURB1,  an  approximation  to  r0  is  obtained  by  integrating 
the  local  asymptotic  value;  however,  with  the  Shkarofsky  model  it  is  a  simple  matter 
to  compute  the  exact  coherence  function. 

IV.  1  Plane  and  Spherical  Waves 

For  a  plane  wave,  the  mutual  coherence  function  is  position-invariant  and  has  already 
been  given  by  (21).  For  a  spherical  wave  emanating  from  the  origin  of  the  reference 
coordinate  system,  it  can  be  shown  that 

T(Ap;  pcl  z)  =  \  exp{t -pc  •  A p  -  tf(A/>;  z)},  (34) 

2*  Z 

where 

H(Ap;z)  =  j  J‘ D,(yz’ I z)dz‘.  (35) 

There  is  a  deterministic  dependence  on  position  relative  to  the  source  that  can  be  impor¬ 
tant.  The  spreading  loss  is  given  by  the  z~2  dependence,  but  this  factor  is  conveniently 
absorbed  in  the  gain  factors  for  the  optical  system  of  interest.  The  most  important  tur¬ 
bulence  effect  for  spherical  waves  is  the  variation  in  the  effective  coherence  scale  along 
the  integration  path — the  shower-glass  effect.  Regions  of  very  high  turbulence  near  the 
source  are  not  as  troublesome  as  regions  of  moderate  turbulence  further  along  the  path. 

IV. 2  Beam  Waves 

The  fundamental  quantity  of  interest  is  the  complex  wavefield  u(r),  but  for  narrow-angle 
scatter  it  is  convenient  to  introduce  the  modified  wavefield 

U(p,  z)  =  -u(r)exp{—  ikz}.  (36) 
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For  a  beam  wave, 


(37) 


U(p,0)  =  «cp{-ifcaf>*}, 
which  has  the  two-dimensional  Fourier  transform 

t7(K,0)  =  i«p(-|l}.  (38) 

Beyond  the  z  =  0  plane, 


U(p>z) 


J f  U(K) 0) exy{ik(g(K)  -  l)z}exp{tK  • 
1  (  ka  p2 

TTi^'xp{-Yrn^}' 


(39) 


To  evaluate  the  integral,  the  narrow-angle  scatter  approximation  g(K)  =  ^Jl  —  (K/k)*  & 
1  —  (K/k)2/ 2  was  used.  It  follows  that 


Mp,  *)!’  =  {Wo/W{z))2  exp{-2  p3/W2(z)},  (40) 


where 

=  +  («.*)*)•  (41) 

Ishimaru  sets  ar  =  2/(&W£)  and  a,*  =  l/J?o,  whereby  W(z)  achieves  a  minimum  at 
z  «  ±1  /Rq.  Thus,  the  beam  wavefield,  which  has  initial  size  W0)  either  converges  to 
a  minimum  at  z  —  1/Ro  (focused)  or  diverges  from  a  fictitious  point  at  z  —  —l/Ro 
(unfocused). 

Under  the  conditions  of  narrow-angle  scatter  and  small  local  perturbations,  the  MCF 
of  an  arbitrary  wavefield  admits  the  integral  representation 

r(Ap,pc,z)  =  y/r.(Ap-zKJ/t;KJ;0)exp{.-KJpe} 

xtxrt-m&p-zKi/k-.KjflJKj,  (42) 

where 

r o{pd,  Kj;  0)  =  J j  T0{pd,  pe\  0)  expl-iKrf  •  pc}  dpc  (43) 

and 

H(Ap;  Kd)=-jJ‘D„(Ap  +  z'Kd/k)  dz'.  (44) 

The  quantity  ro(Ap,  pc,0)  represents  the  MCF  of  the  incident  wavefield  at  z  —  0  in 
terms  of  the  sum  variable  pc  =  (pt  +  pa)/2  and  the  difference  variable  Ap  =  p2  —  p{. 
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For  a  plane  wave,  T0(pj,  pc,  0)  =  1,  which  implies  that 

rO(PdiKd,0)  =  £(Kd) 

independent  of  pd.  Substituting  (45)  into  (42)  gives 


r(Ap,pc,z)  =  e*p{—y  jf  Dn(Ap)dz'}, 


which  is  independent  of  pc.  The  Ap  dependence  of  the  final  result  comes  from  //(A/*;  0). 
For  a  spherical  wave, 


so  that 


It  is  readily  shown  that 


Uo{p,z)  =  -  exp  {ikp2l(‘2z)}, 
z 


r0(Ap,  pc ;  z)  =  -  cxp{ikpc  ■  A p/z). 


T(pJ-,KJ,z)  =  -S(KJ-kpJ/z). 


Taking  the  limit  as  z  — »  0,  it  follows  that 

r(p;K»;0)  =  JjS(pj). 

Substituting  (50)  into  (42)  as  before  gives 

r./A  \  expfi&Ap  •  pjz}  .  h 2 

T(Ap,  pc;z)  =  -AA - f—tsLJ.  exp{  — 


P{~Y  jf  Dn{&pz'/z)dz'}. 


The  leading  term  in  (51)  is  the  MCF  for  the  spherical  wave  source. 

For  a  beam  wave,  it  follows  from  (37)  that 

r0(Ap,  pc,  0)  =  exp{-fc[ar(/)^  +  A p2jA)  +  ia,pc  •  Ap)}. 
Substituting  (52)  into  (43)  it  follows  by  direct  computation  that 


I 'o{Pd\  Kj;  0)  =  exp{—karpd/A  -  [Kj  +  ka{pd}2 /(Akar)}. 


Furthermore, 


lim  r0(pd;  Kd]  0)  =  6(Kd  +  fca,pd). 


If  or,-  ~  0,  we  recover  the  plane- wave  result.  At  the  other  extreme  of  large  a;  we  have 

lim  l'0(pd ;  Kj;  0)  «  j7^~z;S{pd).  (55 

°r  >o  (A 
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IV.3  Small  ar  Approximation 


For  the  moment,  let  us  assume  that  we  are  free  to  make  ar  as  small  as  we  like.  In  that 
case  the  integral  of  (42)  with  T0  defined  by  (53)  can  be  evaluated  by  the  saddle  point 
method.  After  some  straightforward  manipulations  we  find  that 

Y(&p)Pc]z)  «  (W0/W(z))iex?{b2/(4a)-cApi} 

xexp{-yjf  Dn(Ap  +  (z'  -  z^ytydz'},  (56) 


where 


b  =  (2fcfa,(l  -  a;z)  -  a)a)Ap  +  iika.p^  «  a‘*^  Ap,  (58) 

+  “  S'  (59) 


K;  =  — b/(2a) 


A:a,Ap 
1  —  a  ;z 


Upon  substituting  the  approximate  forms  of  (57),  (58),  (59),  and  (60)  into  (56),  the 
dependence  on  a,  cancels  and  the  result  implied  by  (54)  is  obtained.  The  accuracy  of 
the  saddle  point  method  depends  on 

F( Kj)  =  D(Ap  +  w  -  *)( KJ  +  K j)/k)dz'}  (61) 

varying  slowly  over  Kj  regions  such  that  |Kj|  «  1/a.  A  test  of  this  condition  can  be 
made  upon  evaluating  the  integral  in  (56).  The  three  MCF  forms  are  summarized  in 
Table  2. 
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Table  2.  Summary  of  mutual  coherence  function  forms. 


Wave  Type 

r(A  p\pc\z) 

Plane 

exp{-y  Iq  Dn(Ap)  dz '} 

Spherical 

exp{_j£  jj  Dn{Apz'/z)  d2.} 

Beam 

(Wo/PV(z))2  exp{62/(4a)  —  cAp2} 

X  exp{-f  £>„(Ap  +  (z‘  -  z)K’d/k)  dz'} 

IV.4  Intensity  Statistics 

The  fourth-order  moment  of  the  complex  field  satisfies  a  known  differentia)  equation, 
which  has  received  considerable  attention  over  the  past  five  years.  Through  numerical 
simulations,  diagram  methods,  asymptotic  expansions,  and  other  techniques  a  detailed 
understanding  has  emerged,  but  the  most  general  results  cannot  yet  be  encoded  into 
simple  formulas.  For  our  purposes  here,  we  shall  restrict  ourselves  to  single-frequency 
two-point  correlations,  which  are  the  direct  extension  of  the  mutual  coherence  function 
analysis.  For  a  single  pi '  e  screen,  it  can  be  shown  [8,  9]  that  the  spatial  wavenumber 
spectrum  of  the  intensity  k  given  by  the  integral 

i,(K)  =JJ  «p{-S(i,  K! $)}  cos(K  •  0  di  (62) 

where 

=  8k*lp  [J  *(q,0)sin!(i)  ■  q/2)siu!((  •  (63) 

l)  =  z/k,  (64) 

and  lp  is  the  path  length  as  in  Section  III.  The  two-dimensional  Fourier  transform  of 
(63)  gives  the  intensity  correlation  function;  however,  we  shall  be  mainly  interested  in 
the  intensity  moment 

<7W/*‘<K>(f^  <65' 


18 


To  remove  the  mean  intensity,  we  define 


a)  =</2>  -1,  (66) 

which  is  strictly  valid  for  plane  waves  but  a  good  approximation  for  most  applications. 
If  we  substitute  (22)  into  (63)  and  change  variables,  it  can  be  shown  that 

sU.Kl})  =  81 jJJ *in!(Kl,  •  q/2)Sin!(f  •  q/(2i,))  (67) 

where 

1'  =  [k%{Q.0Z3C2n)f3 

V  =  ((///')6/3 

'll  =  Aj// 

'1m  =  (^m^/) 

The  main  point  here  is  that  the  behavior  of  the  intensity  statistics  depends  on  two 
length  parameters,  If,  which  is  the  Fresnel  radius,  and  which  is  close  to  the  receiver 
coherence  diameter,  and  two  ratio  parameters,  7*  and  -ym .  If  the  parameters  7 j  and/or 
7m  become  large,  the  influence  of  the  corresponding  scale  size  becomes  small.  The 
intensity  statistics  are  well  defined  in  the  fractal  limit  where  both  the  outer  and  inner 
scale  sizes  are  respectively  infinity  and  zero.  For  optical  turbulence  effects,  however, 
the  inner  scale  generally  cannot  be  neglected. 

For  the  moment,  let  us  consider  the  fractal  limit.  If  we  replace  /'  with  lc  such  that 
k2lpDn(lc)  as  1,  it  can  be  shown  [8]  that 

where 

U  =  \k%Dn{lf).  (69) 

The  approximation  is  valid  as  long  as  U  «  1.  As  U  increases  beyond  this  range,  cr)  will 
increase  monotonically  until  it  slightly  exceeds  unity  and  then  will  decrease  to  unity, 
its  saturation  value.  The  effects  of  the  inner  scale  are  shown  in  Figure  1  of  Whitman 
and  Beran  [10].  The  general  effect  is  to  shift  the  7m  =  0  curve  to  the  right  and  increase 
its  amplitude.  Simulations  by  Martin  and  Flatte  [11]  have  shown  that  Whitman  and 
Beran’s  calculations  are  about  15%  too  low,  but  the  general  shape  is  preserved.  For 
modeling  purposes,  it  would  not  be  too  difficult  to  fit  a  simple  functional  form  to  the 
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curves  derived  numerically  or  from  simulations.  For  strong  intensity  fluctuation,  we  can 
use  the  saturation  result 

<//'>  -1  «  exp{-2tf},  (70) 

that  is,  we  simply  square  the  mutual  coherence  function.  For  weak  intensity  fluctuations 
approximations  based  on 


*/(K)  *  4 k\  2/}/2)$n(K))  (71) 

can  be  used. 
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V  Slow  and  Fast  Integration  Times 


In  a  turbulent  flow  field,  structures  of  different  sizes  travel  at  different  speeds.  Nonethe¬ 
less,  the  dominant  propagation  effects  come  from  larger  scale  structures  that  are  essen¬ 
tially  frozen  over  the  measurement  integral.  Thus,  a  long  exposure  image  in  the  focal 
plane  of  an  optical  system  will  tend  to  measure  the  ensemble  average  of  the  angular 
deviation  predicted  by  the  MCF.  Shorter  sequential  exposures  tend  to  show  sharper  de¬ 
tail  that  rearrange  themselves  slowly  in  time.  Assuming  that  these  displacements  come 
from  the  larger  scale  turbulent  structures,  one  is  motivated  to  define  the  short-term 
MCF  as  a  high-pass  filtered  version  of  its  long-term  form.  The  concept  is  an  old  one, 
although  it  has  no  firm  theoretical  basis.  The  particular  measure  we  have  implemented 
was  proposed  by  Fante  [12]. 

For  isotropic  turbulence,  the  correlation  function  and  spectral  density  functions  are 
related  by  Hankel  transforms.  Thus,  if 


C(Ap)  =  2x  [°°  KJoiKApfai^dK, 
Jo 

UK)  =  ~  jH  ApC(Ap)J0(KAp)dAp. 


It  follows  that 


C.(Ap)  =  2x  l°°  KJ0(KAp)it(K)dK 

JKq 

=  C(Ap)  -  2?r  /  *  KJ0(KAp)$t(K)dK.  (74) 

Jo 

Because  we  have  an  analytic  form  for  C(Ap ),  it  is  desirable  to  perform  the  filter  oper¬ 
ation  in  the  spatial  domain.  It  can  be  shown  from  (74)  that 

C.(Ap)  =  G{ Ap)  -  [°°  Ap'C(Ap')F(Ap,  Ap')  dAP)  (75) 

Jo 


where 


P(A,n  A  JV  _  tf  r  +  V*  -  ZApAp'cosB)  dQ 

f(AP)Ap)-K0Jo 


'  r  r/  Jo  VAp2  + A/)'l-2A/9Ap/cos0  2ir  v  7 

The  short-term  structure  function  can  then  be  written  as 

D,(AP)  =  Dn(Ap)  -  0.5(C.(0)  -  C,(Ap)).  (77) 

The  short-term  correction  to  be  applied  to  Dn(Ap)  to  obtain  D,(Ap )  is  obtained  by 
interpolating  a  precalculated  set  of  offset  curves.  Prior  to  computing  the  offset  curves, 
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the  structure  functions  are  normalized  to  £7*,  which  effectively  reduces  the  dependence 
on  the  structure  constant  to  a  scale  factor.  It  was  found  that  the  dependence  on  the 
inner  scale  is  very  weak,  so  that  a  nominal  value  of  .5  cm  could  be  used.  The  remaining 
parameters  are  A p)  Lq)  and  which  is  expressed  as  a  fraction  of  2irjLo>  The  effect 
of  changing  Kq  is  similar  to  changing  Lq  itself.  Figure  4  shows  a  family  of  five  D% 
curves  for  Kq  =  0.1,  0.25,  0.5,  2.0,  and  3.0  times  2k/Lq.  As  the  fast  cutoff  wavenumber 
increases,  the  correlation  at  a  fixed  separation  increases.  This  implies  a  sharper  image, 
which  may  be  disp’aced. 
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YI  Simulation  of  Imaging  in  Turbulence 

VI.  1  Background 

The  IMTURB  code  provides  quantitative  system  performance  guidelines  for  imaging 
systems  operating  in  the  turbulent  atmospheric  boundary  layer,  but  it  is  very  difficult  to 
design  definitive  validation  tests  for  the  model.  Ideally,  one  would  design  experiments 
that  combined  channel  diagnostics  with  various  measurements  of  image  degradation. 
Alternatively,  numerical  simulations  can  be  used  to  provide  preliminary  checks  for  uni¬ 
form  propagation  environments.  This  section  describes  a  series  of  measurements  that 
were  performed  to  test  the  short-term  coherence  function  and  beam-wander  concepts. 


VI.2  Propagation  Model 


The  propagation  of  a  spherical  wave  can  be  derived  from  the  relation 


exp{tkr} 
4  7IT 


/Ytexp{ikg(if)lzl} 
JJ  2kg(K) 


exp{iK  •  p}  dK/(27r)2, 


(78) 


which  shows  that  i  exp{ikg(K)\z\} /(2kg(K))  is  the  two-dimensional  Fourier  transform 
of  the  outward-propagating  spherical  wave  in  a  constant  z  plane.  Forward  propagation 
to  any  other  plane  is  achieved  by  applying  the  propagation  factor  exp{ikg(K)Az}, 
which  simply  extends  the  z  component.  Any  attempt  to  approximate  this  relation 
numerically,  however,  is  futile  because  of  the  slow  decay  of  the  field.  The  integral  itseH 
is  defined  only  in  a  special  mathematical  sense. 

The  problem  can  be  avoided  by  using  a  divergent  beam,  which  retains  the  phase 
variation  spherical  wave  but  imposes  a  gaussian  intensity  profile.  For  a  beam  wave, 


Mp>z) 


i 

1  +  iaz 


exp {ikz  — 


2  l  +  iaz*' 


where  a  =  aT  +  ia,-  with 


(79) 


a,  = 

a;  = 


(80) 

(81) 


As  discussed  previously,  the  beam  size  at  z  =  0  is  W0>  and  Rq  is  the  curvature  of  the 
phase  front.  For  a  focused  beam,  the  beam  intensity  spread  reaches  a  minimum  size  u>0 
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at  z  «  Re-  The  beam  size  as  a  function  of  distance  from  this  minimum  is  given  by  the 
relation  _ 

W(z)  =  u0Jl  +  -^~*z*.  (82) 

V  *  wo 

At  large  distances,  the  angular  extent  of  the  beam  is  \/(nu>0).  Thus,  the  beam  is  well 
approximated  numerically  if  the  spatial  wavenumber  at  the  maximum  angular  deviation 
is  resolved.  This  leads  to  the  sampling  condition 

Ax  <<  ttu>o/2.  (83) 

The  number  of  points  needed  is  determined  by  the  size  of  the  beam  at  the  position 
where  the  propagation  disturbance  is  applied. 

Following  Martin  and  Flatte  [4]  we  simulate  the  effects  of  strong  turbulence  by  multi¬ 
plying  the  beam  field  by  a  complex  random  field  whose  phase  autocorrelation  function 
is  given  by  12  and  13  using  the  model  summarized  in  Table  1.  The  location  of  the 
phase  screen  is  chosen  so  that  the  beam  is  well  tapered  over  the  finite  extent  of  the 
phase  screen.  We  then  propagate  the  distorted  beam  wavefield  forward  by  performing 
a  Fourier  transform  and  then  multiplying  each  Fourier  component  by  the  propagation 
factor  exp{ikg(K )Az}.  Strict  adherence  to  the  parabolic  wave  equation  requires  repeti¬ 
tion  of  this  process,  but  experience  has  shown  that  the  field  statistics  well  removed  from 
the  phase  screen  are  not  significantly  different  if  a  single  phase  screen  corresponding  to 
the  total  path-integrated  disturbance  is  used. 

In  Martin  and  Flatte’s  work,  multiple  phase  screens  were  used,  but  they  simulated 
only  the  intensity  of  the  scattered  field.  To  simulate  the  blurring  effects  of  the  turbulence 
in  an  image,  we  have  numerically  refocused  the  beam  to  z.  This  is  accomplished  by 
back  propagating  the  disturbed  beam  wavefield  through  the  phase  screen  location  to 
the  point  of  minimum  beam  diameter.  The  same  FFT  procedure  is  used  to  implement 
the  back  propagation.  In  the  absence  of  a  phase  disturbance,  we  obtain  a  small  spot  as 
would  be  observed  with  a  telescope  pointed  along  a  focused  laser  beam.  Thus,  we  use 
a  focused  1.064  /im-beam  with  a  1-cm  waist  diverging  to  approximately  30  cm  at  the 
phase  screen.  The  phase  screen  was  placed  approximately  2/3  of  the  distance  from  the 
beam  waist  to  the  aperture  plane. 

Figures  5,  6,  and  7  show  the  effects  of  weak,  moderate,  and  severe  turbulence, 
respectively.  The  structure  constant  and  the  path  length  within  the  disturbed  medium 
are  listed  in  the  figures.  The  distances  are  real  units  at  the  location  of  the  beam 
minimum.  Each  realization  of  the  focused  image  can  be  thought  of  as  an  instantaneous 
field  sample  as  could  be  obtained,  for  example,  by  a  video  tape  recording  of  the  beam 
focused  on  a  screen.  Long-term  averages  are  obtained  by  incoherently  averaging  a  large 
number  of  realizations.  Figures  1,  2,  and  3  suggest,  however,  that  the  beam  centroid  is 
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meters 


Figure  5.  Focused  beam  simulation  for  weak  turbulence. 


near  tlie  origin,  with  very  little  deviation.  We  shall  6ee  in  the  next  section  that  this 
indeed  the  case. 
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Figure  6.  Focused  beam  simulation  for  moderate  turbulence. 
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Gprop  Run  06  Cn**2=1.0e-12  nx=ny=64  lo=0.005m  Lo=1m  Lp= 1000m 


X  -  meters  Wmin=0.01m  Wscr=nx*dx/4  W(z=0)=10*resolution 


Figure  7.  Focused  beam  simulation  for  strong  turbulence. 
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VI.3  Summary  Statistics 

To  summarize  the  average  blurring  effects  or  long-time  images,  we  generated  a  large 
number  of  independent  realizations  for  several  different  turbulence  levels  and  computed 
standard  measures  of  image  degradation.  The  first  moment  and  second  central  mo¬ 
ment  of  the  image  intensity  along  the  reference  coordinate  directions  measure  average 
displacement  and  image  spread,  respectively.  We  also  computed  the  autocorrelation 
function  of  each  image  as  a  measure  of  its  fine  structure.  Finally,  we  computed  the 
autocorrelation  function  of  the  average  intensity  of  a  large  number  of  independent  re¬ 
alizations  as  a  measure  of  the  long-term  average  spot  size. 

As  a  measure  of  the  potential  effects  of  the  turbulence,  we  compute  the  reciprocal  of 
the  e-1  decorrelation  distance  at  the  aperture  plane  for  the  theoretical  mutual  coherence 
function.  In  the  IMTURB  code,  the  coherent  spot  diameter  is  inversely  proportional 
to  the  coherence  scale.  Thus,  we  have  presented  the  simulation  results  in  terms  of  the 
reciprocal  coherence  scale.  To  provide  a  unitless  measure,  the  coherence  scale  should 
be  normalized  to  the  aperture  size;  however,  only  a  single  aperture  size  was  used  for 
this  set  of  simulations. 

Figures  8,  9,  10,  and  11  show  the  average  of  the  measured  first  and  second  moments. 
Because  the  turbulence  is  isotropic,  we  expect  and  observe  no  average  displacement  or 
difference  between  the  *  and  y  cuts.  Although  the  spread  in  the  displacement  increases 
with  increasing  perturbation  strength,  the  wander  is  much  smaller  than  the  average 
spread.  Thus,  if  the  images  were  displayed  sequentially,  no  beam  wander  would  be 
apparent.  In  rl,ort,  the  homogeneous  turbulence  model  does  not  contain  sufficient  low 
spatial-frequency  content  to  account  for  the  beam  wander  that  is  typically  observed  in 
experiments. 

The  most  interesting  result  is  summarized  in  Figure  12,  which  shows  the  average 
short-  and  long-term  structure  measures.  The  long-term  beam  size  starts  at  twice  the 
focused  spot  size  because  of  the  autocorrelation  function  measure.  As  the  turbulence 
develops,  however,  the  long-term  beam  size  increases  linearly  against  the  inverse  coher¬ 
ence  scale,  as  do  the  second  moments  shown  in  Figures  10  and  11.  This  confirms  the 
known  inverse  proportionality  between  the  average  focused  spot  size  and  the  coherence 
scale  at  the  aperture  plane. 

The  increase  in  the  short-term  structure  size,  however,  is  not  linear  with  inverse 
coherence  scale.  This  is  possibly  due  to  Fresnel-radius  effects,  which  are  intimately  part 
of  intensity  fluctuations.  In  any  case,  we  see  that  the  short-term  image  coherence  scale 
increases  at  an  increasing  rate  as  a  function  of  inverse  coherence  scale.  It  is  well  known 
that  short-term  exposures  can  be  refocused  to  remove  the  blurring  with  adaptive  phase 
compensation  schemes.  We  suspect  that  the  short-term  coherence  scale  presented  here 
represents  the  best  that  can  be  achieved  with  such  schemes.  For  example,  in  very  strong 
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Beam  Propagation  Effects  wavelen=  1.064  micron  50  surf 


1/lc  ap.  -  meier-l  wo=0.14!42m  wmin=0.01m  Lo=lm  !o=.005m 

Figure  8.  Average  beam  displacement  along  x-axis. 

turbulence,  adaptive  focusing  can  recover  only  about  50%  of  the  diffraction  limited 
resolution.  If  this  can  be  confirmed  by  simulating  adaptive  phase  reconstruction,  it  will 
provide  an  improved  quantitative  measure  of  the  short-term  spot  size  currently  provided 
by  the  IMTURB  code. 
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Figure  9.  Average  beam  displacement  along  y-axis. 
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Figure  12.  Long-  and  short-time  image  correlation  measures. 
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VII  Conclusions  and  Recommendations 


A  consistently  formulated  propagation  model  has  been  described  for  systems-oriented 
predictive  codes.  The  primary  model  inputs — the  structure  constant  and  the  outer  and 
inner  scale  sizes — can  be  obtained  from  atmospheric  turbulence  models.  In  terms  of 
simple  path  integrals,  the  complete  mutual  coherence  function  for  plane-waves,  point 
sources,  and  beams  can  be  computed.  Iu  fact,  the  plane  and  point  source  limits  are 
special  cases  of  the  beam  wave  result.  The  model  eliminates  the  restrictions  imposed 
by  the  Rytov  approximation. 

Because  of  the  importance  of  image  processing  and  adaptive  compensation  of  at¬ 
mospheric  turbulence,  we  also  implemented  a  consistent  measure  of  the  short-term 
coherence  scale.  As  with  all  previous  measures,  however,  there  is  no  firm  theoretical 
foundation  for  the  result.  Thus,  we  used  numerical  simulations  to  study  the  relation 
between  the  long-term  and  short-term  correlation  scales.  Moreover,  in  our  introduc¬ 
tory  background  material,  we  showed  that  any  model  based  on  the  assumption  of  locally 
homogeneous  statistics  eliminates  phase  wander  that  is  inherently  part  of  any  real  prop¬ 
agation  environment.  Our  simulations  confirm  that  slow  beam  wander  is  not  part  of 
standard  homogeneous  turbulence  models.  To  the  extent  that  it  needs  to  be  included 
in  predictive  models,  it  must  be  accommodated  explicitly. 

More  importantly,  we  developed  a  quantitative  measure  of  the  short-term  intensity 
coherence.  We  don’t  know  how  to  calculate  the  measure  theoretically,  but  it  could 
easily  be  modeled  empirically.  As  noted  in  Section  VI,  we  believe  that  the  short-term 
coherence  scale  as  we  defined  it  is  a  quantitative  measure  of  what  can  be  achieved 
by  adaptive  focusing.  Future  work  should  test  this  concept.  Beyond  that,  further 
improvements  to  the  IMTURB  code  should  be  coordinated  with  field  measurements. 
In  a  separate  memorandum,  we  have  described  a  simple  and  cost-effective  laser  diagnostic 
system  that  could  achieve  this  objective. 
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Appendix 


Optical  Systems 

Ultimately,  the  wave  fields  that  have  been  characterized  in  this  report  are  detected  and 
processed  by  an  optical  system.  The  performance  of  the  optical  system  for  its  intended 
purpose  is  critically  dependent  on  its  detailed  configuration,  but  there  are  some  aspects 
that  are  common  to  all  systems.  In  this  appendix  a  general  model  that  isolates  these 
common  features  is  described. 


A-I  A  Model  for  Optical/IR  Imaging 

Consider  the  general  model  shown  in  Figure  Al.  Any  luminous  object  can  be  decom¬ 
posed  into  a  collection  of  effective  point  radiators,  but  for  simplicity  let  0(£)  represent 
a  collection  of  noninteracting  point  radiators  that  are  visible  to  the  optical  system.  Let 
A Sr(p,£)  represent  the  complex  wave  field  in  the  aperture  plane  of  the  receiver  from  a 
unit  point  radiator  at  £.  For  example,  in  free  space, 

ASS(,,0  =  SS^,  (i.) 

T 

where 

r  =  jB?  +  2R-(p-()  -l-lp-il2-  (42) 

It  follows  from  the  linearity  of  Maxwell’s  equations  that  the  aperture-plane  signal  from 
any  extended  luminous  object  can  be  represented  by  the  integral 

Sr(j>)=  I,  ASR(p,()0(()d(.  (43) 

•/object 

If  the  luminosity  is  broadband,  (.43)  must  be  integrated  over  the  temporal  frequency 
spectrum,  but  the  narrowband  model  is  adequate  to  illustrate  the  linear  systems  ap¬ 
proach.  The  complex  signal  Sr(p)  over  some  finite  aperture  A  is  accessible  to  the 
imaging  system  for  processing.  The  form  of  (A3)  is  the  general  representation  of  a 
spatially  varying  linear  system. 

The  simplest  optical  system  uses  a  collection  of  lenses  that  effectively  apply  phase 
compensation  for  the  propagation  paths  to  p  from  each  point  in  the  object  space,  £o- 
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Figure  Al.  Geometry  for  general  imaging  model. 


This  operation  can  be  represented  mathematically  as 

=  /  SR(p)  exp{t^H(/>,  &)}  dp 

J  A 

=  /  /  ^-?H(p,0^(0^exP{iV’/i(/,.^o)}dp.  (Al) 

J  A  J  object 


The  second  form  is  obtained  by  substitution  from  (v43).  Changing  the  order  of  integra¬ 
tion  in  (i44)  gives  the  linear  systems  representation 


V«)=/ 

J  object 

where 

M(t,io)=  [  &SR(p)t)exi>{i'ij>R(p,to)}dp 
J  A 


(45) 

(46) 
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is  the  complex  point  spread  function  for  the  imaging  system.  For  an  ideal  system, 
M((,h)  —  £(£  -  £o)  times  an  arbitrary  phase  constant. 

Over  at  least  a  limited  portion  of  the  receiving  aperture,  the  point  spread  function  of 
a  good  imaging  system  will  depend  only  on  the  difference  between  the  object  and  target 
coordinates.  Thus,  the  first  step  in  any  imaging  process  is  mathematically  equivalent 
to  applying  a  homogeneous  spatial  filter  to  the  object  illumination.  The  differences 
between  coherent  and  incoherent  systems  manifest  themselves  in  the  characteristics  of 
the  image  intensity,  particularly  its  temporal  and  spatial  coherence  properties.  In  all 
cases,  however,  if  propagation  disturbance  are  present,  the  point  response  function 
itself  acquires  a  random  time- varying  component. 


A-II  Propagation  Disturbances  in  Optical  Systems 

To  identify  the  random  component  induced  by  the  propagation  disturbances  explicitly, 
let 

ASr(p,()~Hk(p,()ASUp,().  (ai) 

Now  consider  an  object  whose  luminosity  is  uncorrelated  from  point  to  point  and  inde¬ 
pendent  of  the  propagation  disturbance.  It  follows  from  (45)  that 

C|V(WI!>=  /  <|M«,Wr><|0(0|!>  i(,  (AS) 

•/object 

where  the  angle  brackets  denote  ensemble  averaging.  Similarly,  from  (46)  it  follows 
that 

=  JJa  <Mp,omp‘,o> 

x  exp £o)  ~  ^r(p',  &>)]}  dpdp'.  (49) 

Because  R  »  p  or  £,  (42)  can  be  expanded  to  terms  that  are  quadratic  in  R~l, 
whereby 

&sr{p>  0&sr*(p'> 0  ~  exp{-ikAp  ■  i/R},  (410) 

to  terms  that  do  not  depend  on  £  or  are  quadratic  in  £/\/R.  For  an  imaging  system, 
there  is  a  similar  relation  for  ipR  with  R  replaced  by  2/,  where  /  is  the  focal  length. 
Also,  in  Section  III.  1  we  showed  that 

r(A/>;{,  fi)  =<f Up, Ome'.O  a s%(p,()as°b-(p',(),  (a n) 

so  that 

<Hr(p>0H'r(p',0>=  exp{-//(A  />,£; R)}.  (412) 
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Here  we  have  allowed  for  a  weak  dependence  of  H(Ap,  R)  on  If  this  dependence 
is  neglected,  the  propagation  disturbance  is  isoplanatic.  Finally,  we  define  the  function 
W{p )  to  be  unity  on  A  and  zero  elsewhere  and  introduce  the  new  variables 

v  =  H/R 

Vo  =  *£o/(2/). 

With  these  definitions  and  approximations,  it  follows  that 

<M(t],t] 0)>«  jF(Ap)exp{-H(Ap,t]R/k’,R)}exi>{-iAp-Aj]}dAp,  (>415) 

where 

F(Ap)  =  J\V(x+  A p)W (x  -  A p)  dx-  (>116) 

If  there  is  no  propagation  disturbance,  H  =  1,  and  from  (>115)  we  see  that  the  point 
spread  function  M(Ar])}  which  is  no  longer  random,  and  F(Ap)  form  a  Fourier  transform 
pair.  Analogous  to  the  impulse  response-transfer  function  relation  for  linear  filters, 
F(Ap)  is  called  the  modulation  transfer  function  (MTF)  for  the  optical  system.  Note 
that  Goodman  [13]  uses  the  terminology  optical  transfer  function.  If  the  disturbance  is 
isoplanatic,  we  can  define  an  average  MTF  as 

<F(Ap)>-  F{Ap)exp{-H{Ap]R)}.  (A17) 

Equation  (A17)  is  the  long-exposure  MTF  for  an  optical  system — see  Fried’s  equations 
(3.15)  and  (3.16)[14]. 

With  short  exposures  it  is  well  known  that  one  obtains  a  sharper  image  that  is 
displaced  from  its  center  point.  It  is  mainly  the  slow  meandering  of  the  diffracted  beam 
that  contributes  to  the  blurring  of  a  point  image.  To  estimate  the  short  exposure  MTF, 
it  is  clear  that  one  should  apply  a  high-pass  spatial  filter  to  (A 7)  before  performing  the 
ensemble  average  in  (A6).  Unfortunately,  this  does  not  yield  mathematically  tractable 
results,  and  various  other  schemes  have  been  used.  Fante  [15,  12],  for  example,  has 
modified  the  lower  limit  of  int  egration  in  (17)  from  0  to  -yAp/D ,  where  7  «  1,  and  D  is 
the  size  of  the  largest  unresolved  eddy,  but  he  remains  skeptical  of  the  procedure  [16]. 
Fried  [14]  removes  the  least-squares  estimate  of  the  linear-phase  component;  however, 
the  general  area  of  tilt  correction  merits  a  more  thorough  treatment.  At  the  present 
time,  however,  there  is  no  satisfactory  analytical  approximation  for  the  short-term  MTF. 


(413) 

(414) 
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